#pragma once



class KisFourPointInterpolatorBackward
{
public:
    KisFourPointInterpolatorBackward(const QPolygonF& srcPolygon, const QPolygonF& dstPolygon) {
        m_a = dstPolygon[1] - dstPolygon[0]; // AB
        m_b = dstPolygon[2] - dstPolygon[1]; // BD
        m_c = dstPolygon[3] - dstPolygon[0]; // AC
        m_d = m_b - m_c; // BD - AC

        m_qA = m_c.x() * m_d.y() - m_c.y() * m_d.x();

        m_srcBase = srcPolygon[0];
        m_dstBase = dstPolygon[0];
        m_xCoeff = srcPolygon[1].x() - srcPolygon[0].x(); // AB_src
        m_yCoeff = srcPolygon[3].y() - srcPolygon[0].y(); // AC_src

        m_qB_const = m_c.x() * m_a.y() - m_c.y() * m_a.x();

        m_qD_div = 1.0 / (2 * m_qA);

        //m_qB_varX = 0.0;
        //m_qB_varY = 0.0;
    }

    inline QPointF map(const QPointF& pt) {
        setX(pt.x());
        setY(pt.y());
        return getValue();
    }

    inline void setX(qreal x) {
        x -= m_dstBase.x();

        m_qB_varX = -x * m_d.y();
        m_qC_varX = -x * m_a.y();
        m_px = x;
    }

    inline void setY(qreal y) {
        y -= m_dstBase.y();

        m_qB_varY = y * m_d.x();
        m_qC_varY = y * m_a.x();
        m_py = y;
    }

    inline QPointF getValue() const {
        static const qreal eps = 1e-10;

        qreal qB = m_qB_const + m_qB_varX + m_qB_varY;
        qreal qC = m_qC_varX + m_qC_varY;

        qreal nu = 0.0;

        if (qAbs(m_qA) < eps) {
            nu = -qC / qB;
        }
        else {
            qreal D = qB * qB - 4 * m_qA * qC;
            if (D > 0.0) {
                qreal sqrtD = std::sqrt(D);
                nu = (-qB - sqrtD) * m_qD_div;
                if (nu < 0.0 || nu > 1.0) {
                    qreal nu2 = (-qB + sqrtD) * m_qD_div;

                    if (nu2 < 0.0 || nu2 > 1.0) {
                        nu = qBound(qreal(0.0), nu, qreal(1.0));
                    }
                    else {
                        nu = nu2;
                    }
                }
            }
            else {
                nu = 0.0;
            }
        }

        qreal xBasedDenominator = m_a.x() + nu * m_d.x();

        qreal mu;

        if (qAbs(xBasedDenominator) > eps) {
            mu = (m_px - nu * m_c.x()) / xBasedDenominator;
        }
        else {
            mu = (m_py - nu * m_c.y()) / (m_a.y() + nu * m_d.y());
        }

        return m_srcBase + QPointF(mu * m_xCoeff, nu * m_yCoeff);
    }

private:
    QPointF m_a; // AB
    QPointF m_b; // BD
    QPointF m_c; // AC
    QPointF m_d; // m_b - m_c

    qreal m_qA{ 0.0 }; // quadratic equation A coeff
    qreal m_qB_const{ 0.0 }; // quadratic equation B coeff, const part
    qreal m_qB_varX{ 0.0 }; // quadratic equation B coeff, X-dep part
    qreal m_qB_varY{ 0.0 }; // quadratic equation B coeff, Y-dep part
    qreal m_qC_varX{ 0.0 }; // quadratic equation C coeff, X-dep part
    qreal m_qC_varY{ 0.0 }; // quadratic equation C coeff, Y-dep part
    qreal m_qD_div{ 0.0 }; // inverted divisor of the quadratic equation solution
    qreal m_px{ 0.0 }; // saved relative X coordinate
    qreal m_py{ 0.0 }; // saved relative Y coordinate

    QPointF m_srcBase;
    QPointF m_dstBase;
    qreal m_xCoeff{ 0.0 };
    qreal m_yCoeff{ 0.0 };
};


namespace GridIterationTools {



    inline void adjustAlignedPolygon(QPolygonF& polygon)
    {
        static const qreal eps = 1e-5;
        static const  QPointF p1(eps, 0.0);
        static const  QPointF p2(eps, eps);
        static const  QPointF p3(0.0, eps);

        polygon[1] += p1;
        polygon[2] += p2;
        polygon[3] += p3;
    }
    inline QVector<int> calculateCellIndexes(int col, int row, const QSize& gridSize)
    {
        const int tl = col + row * gridSize.width();
        const int tr = tl + 1;
        const int bl = tl + gridSize.width();
        const int br = bl + 1;

        QVector<int> cellIndexes;
        cellIndexes << tl;
        cellIndexes << tr;
        cellIndexes << br;
        cellIndexes << bl;

        return cellIndexes;
    }

    template <class PolygonOp, class IndexesOp>
    struct AlwaysCompletePolygonPolicy {

        static inline bool tryProcessPolygon(int col, int row,
            int numExistingPoints,
            PolygonOp& polygonOp,
            IndexesOp& indexesOp,
            const QVector<int>& polygonPoints,
            const QVector<QPointF>& originalPoints,
            const QVector<QPointF>& transformedPoints)
        {
            Q_UNUSED(col);
            Q_UNUSED(row);
            Q_UNUSED(polygonOp);
            Q_UNUSED(indexesOp);
            Q_UNUSED(polygonPoints);
            Q_UNUSED(originalPoints);
            Q_UNUSED(transformedPoints);

         //   KIS_ASSERT_RECOVER_NOOP(numExistingPoints == 4);
            return false;
        }
    };


    template <template <class PolygonOp, class IndexesOp> class IncompletePolygonPolicy,
        class PolygonOp,
        class IndexesOp>
        void iterateThroughGrid(PolygonOp& polygonOp,
            IndexesOp& indexesOp,
            const QSize& gridSize,
            const QVector<QPointF>& originalPoints,
            const QVector<QPointF>& transformedPoints)
    {
        QVector<int> polygonPoints(4);

        for (int row = 0; row < gridSize.height() - 1; row++) {
            for (int col = 0; col < gridSize.width() - 1; col++) {
                int numExistingPoints = 0;

                polygonPoints = indexesOp.calculateMappedIndexes(col, row, &numExistingPoints);

                if (!IncompletePolygonPolicy<PolygonOp, IndexesOp>::
                    tryProcessPolygon(col, row,
                        numExistingPoints,
                        polygonOp,
                        indexesOp,
                        polygonPoints,
                        originalPoints,
                        transformedPoints)) {

                    QPolygonF srcPolygon;
                    QPolygonF dstPolygon;

                    for (int i = 0; i < 4; i++) {
                        const int index = polygonPoints[i];
                        srcPolygon << originalPoints[index];
                        dstPolygon << transformedPoints[index];
                    }

                    adjustAlignedPolygon(srcPolygon);
                    adjustAlignedPolygon(dstPolygon);

                    polygonOp(srcPolygon, dstPolygon);
                }
            }
        }
    }


    struct RegularGridIndexesOp {

        RegularGridIndexesOp(const QSize& gridSize)
            : m_gridSize(gridSize)
        {
        }

        inline QVector<int> calculateMappedIndexes(int col, int row,
            int* numExistingPoints) const {

            *numExistingPoints = 4;
            QVector<int> cellIndexes =
                GridIterationTools::calculateCellIndexes(col, row, m_gridSize);

            return cellIndexes;
        }

        inline int tryGetValidIndex(const QPoint& cellPt) const {
            Q_UNUSED(cellPt);

         //   KIS_ASSERT_RECOVER_NOOP(0 && "Not applicable");
            return -1;
        }

        inline QPointF getSrcPointForce(const QPoint& cellPt) const {
            Q_UNUSED(cellPt);

       //     KIS_ASSERT_RECOVER_NOOP(0 && "Not applicable");
            return QPointF();
        }

        inline const QPolygonF srcCropPolygon() const {
        //    KIS_ASSERT_RECOVER_NOOP(0 && "Not applicable");
            return QPolygonF();
        }

        QSize m_gridSize;
    };


    template <class ProcessPolygon, class ForwardTransform>
    struct CellOp
    {
        CellOp(ProcessPolygon& _polygonOp, ForwardTransform& _transformOp)
            : polygonOp(_polygonOp),
            transformOp(_transformOp)
        {
        }

        inline void processPoint(int col, int row,
            int prevCol, int prevRow,
            int colIndex, int rowIndex) {

            QPointF dstPosF = transformOp(QPointF(col, row));
            currLinePoints << dstPosF;

            if (rowIndex >= 1 && colIndex >= 1) {
                QPolygonF srcPolygon;

                srcPolygon << QPointF(prevCol, prevRow);
                srcPolygon << QPointF(col, prevRow);
                srcPolygon << QPointF(col, row);
                srcPolygon << QPointF(prevCol, row);

                QPolygonF dstPolygon;

                dstPolygon << prevLinePoints.at(colIndex - 1);
                dstPolygon << prevLinePoints.at(colIndex);
                dstPolygon << currLinePoints.at(colIndex);
                dstPolygon << currLinePoints.at(colIndex - 1);

                polygonOp(srcPolygon, dstPolygon);
            }

        }

        inline void nextLine() {
            std::swap(prevLinePoints, currLinePoints);

            // we are erasing elements for not free'ing the occupied
            // memory, which is more efficient since we are going to fill
            // the vector again
            currLinePoints.erase(currLinePoints.begin(), currLinePoints.end());
        }

        QVector<QPointF> prevLinePoints;
        QVector<QPointF> currLinePoints;
        ProcessPolygon& polygonOp;
        ForwardTransform& transformOp;
    };

    template <class ProcessCell>
    void processGrid(ProcessCell& cellOp,
        const QRect& srcBounds,
        const int pixelPrecision)
    {
        if (srcBounds.isEmpty()) return;

        const int alignmentMask = ~(pixelPrecision - 1);

        int prevRow = std::numeric_limits<int>::max();
        int prevCol = std::numeric_limits<int>::max();

        int rowIndex = 0;
        int colIndex = 0;

        for (int row = srcBounds.top(); row <= srcBounds.bottom();) {
            for (int col = srcBounds.left(); col <= srcBounds.right();) {

                cellOp.processPoint(col, row,
                    prevCol, prevRow,
                    colIndex, rowIndex);

                prevCol = col;
                col += pixelPrecision;
                colIndex++;

                if (col > srcBounds.right() &&
                    col <= srcBounds.right() + pixelPrecision - 1) {

                    col = srcBounds.right();
                }
                else {
                    col &= alignmentMask;
                }
            }

            cellOp.nextLine();
            colIndex = 0;

            prevRow = row;
            row += pixelPrecision;
            rowIndex++;

            if (row > srcBounds.bottom() &&
                row <= srcBounds.bottom() + pixelPrecision - 1) {

                row = srcBounds.bottom();
            }
            else {
                row &= alignmentMask;
            }
        }
    }

    template <class ProcessPolygon, class ForwardTransform>
    void processGrid(ProcessPolygon& polygonOp, ForwardTransform& transformOp,
        const QRect& srcBounds, const int pixelPrecision)
    {
        CellOp<ProcessPolygon, ForwardTransform> cellOp(polygonOp, transformOp);
        processGrid(cellOp, srcBounds, pixelPrecision);
    }






    /*************************************************************/
    /*      Iteration through precalculated grid                 */
    /*************************************************************/

    /**
     *    A-----B         The polygons will be in the following order:
     *    |     |
     *    |     |         polygon << A << B << D << C;
     *    C-----D
     */



     /**
      * There is a weird problem in fetching correct bounds of the polygon.
      * If the rightmost (bottommost) point of the polygon is integral, then
      * QRectF() will end exactly on it, but when converting into QRect the last
      * point will not be taken into account. It happens due to the difference
      * between center-point/topleft-point point representation. In many cases
      * the latter is expected, but we don't work with it in Qt/Krita.
      */


    struct QImagePolygonOp
    {
        QImagePolygonOp(const QImage& srcImage, QImage& dstImage,
            const QPointF& srcImageOffset,
            const QPointF& dstImageOffset)
            : m_srcImage(srcImage), m_dstImage(dstImage),
            m_srcImageOffset(srcImageOffset),
            m_dstImageOffset(dstImageOffset),
            m_srcImageRect(m_srcImage.rect()),
            m_dstImageRect(m_dstImage.rect())
        {
        }

        void operator() (const QPolygonF& srcPolygon, const QPolygonF& dstPolygon) {
            this->operator() (srcPolygon, dstPolygon, dstPolygon);
        }

        void operator() (const QPolygonF& srcPolygon, const QPolygonF& dstPolygon, const QPolygonF& clipDstPolygon) {
            QRect boundRect = clipDstPolygon.boundingRect().toAlignedRect();
            KisFourPointInterpolatorBackward interp(srcPolygon, dstPolygon);

            for (int y = boundRect.top(); y <= boundRect.bottom(); y++) {
                interp.setY(y);

                for (int x = boundRect.left(); x <= boundRect.right(); x++) {

                    QPointF srcPoint(x, y);
                    if (clipDstPolygon.containsPoint(srcPoint, Qt::OddEvenFill)) {

                        interp.setX(srcPoint.x());
                        QPointF dstPoint = interp.getValue();

                        // about srcPoint/dstPoint hell please see a
                        // comment in PaintDevicePolygonOp::operator() ()

                        srcPoint -= m_dstImageOffset;
                        dstPoint -= m_srcImageOffset;

                        QPoint srcPointI = srcPoint.toPoint();
                        QPoint dstPointI = dstPoint.toPoint();

                        if (!m_dstImageRect.contains(srcPointI)) continue;
                        if (!m_srcImageRect.contains(dstPointI)) continue;


                        m_dstImage.setPixel(srcPointI, m_srcImage.pixel(dstPointI));
                    }
                }
            }


            //QPainter gc(&m_dstImage);
            //gc.setPen(Qt::red);
            //gc.setOpacity(0.5);

            //gc.setBrush(Qt::green);
            //gc.drawPolygon(clipDstPolygon.translated(-m_dstImageOffset));

            //gc.setBrush(Qt::blue);
            //gc.drawPolygon(dstPolygon.translated(m_dstImageOffset));



        }

        const QImage& m_srcImage;
        QImage& m_dstImage;
        QPointF m_srcImageOffset;
        QPointF m_dstImageOffset;

        QRect m_srcImageRect;
        QRect m_dstImageRect;
    };

}

struct RegularGridIndexesOp {

    RegularGridIndexesOp(const QSize& gridSize)
        : m_gridSize(gridSize)
    {
    }

    inline QVector<int> calculateMappedIndexes(int col, int row,
        int* numExistingPoints) const {

        *numExistingPoints = 4;
        QVector<int> cellIndexes =
            GridIterationTools::calculateCellIndexes(col, row, m_gridSize);

        return cellIndexes;
    }

    inline int tryGetValidIndex(const QPoint& cellPt) const {
        Q_UNUSED(cellPt);

      //  KIS_ASSERT_RECOVER_NOOP(0 && "Not applicable");
        return -1;
    }

    inline QPointF getSrcPointForce(const QPoint& cellPt) const {
        Q_UNUSED(cellPt);

      //  KIS_ASSERT_RECOVER_NOOP(0 && "Not applicable");
        return QPointF();
    }

    inline const QPolygonF srcCropPolygon() const {
     //   KIS_ASSERT_RECOVER_NOOP(0 && "Not applicable");
        return QPolygonF();
    }

    QSize m_gridSize;
};
